Prom Electrons to Finite Elements: A Concurrent Multiscale Approach for Metals 
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We present a multiscale modeling approach that concurrently couples quantum mechanical, classi- 
cal atomistic and continuum mechanics simulations in a unified fashion for metals. This approach is 
particular useful for systems where chemical interactions in a small region can affect the macroscopic 
properties of a material. We discuss how the coupling across different scales can be accomplished 
efficiently, and we apply the method to multiscale simulations of an edge dislocation in aluminum 
in the absence and presence of H impurities. 



Some of the most fascinating problems in all fields of 
science involve multiple spatial and/or temporal scales: 
processes that occur at a certain scale govern the behav- 
ior of the system across several (usually larger) scales. In 
the context of materials science, the ultimate microscopic 
constituents of materials are ions and valence electrons; 
interactions among them at the atomic level determine 
the behavior of the material at the macroscopic scale, the 
latter being the scale of interest for technological applica- 
tions. Conceptually, two categories of multiscale simula- 
tions can be envisioned, sequential, consisting of passing 
information across scales, and concurrent, consisting of 
seamless coupling of scales The majority of multi- 
scale simulations that are currently in use are sequential 
ones, which are effective in systems where the different 
scales are weakly coupled. For systems whose behavior 
at each scale depends strongly on what happens at the 
other scales, concurrent approaches are usually required. 
In contrast to sequential approaches, concurrent simula- 
tions are still relatively new and only a few models have 
been developed to date 0, S H S S @ . 

A successful concurrent multiscale method is the Qua- 
sicontinuum (QC) method originally proposed by Tad- 
mor et al. |2(. The idea underlying this method is that 
atomistic processes of interest often occur in very small 
spatial domains while the vast majority of atoms in the 
material behave according to well-established continuum 
theories. To exploit this fact, the QC method retains 
atomic resolution only where necessary and grades out 
to a continuum finite element description elsewhere. The 
original formulation of QC was limited to classical po- 
tentials for describing interactions between atoms. Since 
many materials properties depend explicitly on the be- 
havior of electrons, such as bond breaking/forming at 
crack tips or defect cores, chemical reactions with im- 
purities, and surface reactions and reconstructions, it is 
desirable to incorporate appropriate quantum mechani- 
cal descriptions into the QC formalism. In this Letter, 
we extend the original QC approach so that it can be 
directly coupled with quantum mechanical calculations 
based on density functional theory (DFT) for metallic 



systems. We refer to the new approach as QCDFT. 

The goal of the QC method is to model an atom- 
istic system without explicitly treating every atom in the 
problem 0, 0- This is achieved by replacing the full 
set of N atoms with a small subset of N r "representa- 
tive atoms" or repatoms (N r <C N) that approximate the 
total energy through appropriate weighting. The ener- 
gies of individual repatoms are computed in two different 
ways depending on the deformation in their immediate 
vicinity. Atoms experiencing large deformation gradients 
on an atomic-scale are computed in the same way as in a 
standard fully-atomistic method. In QC these atoms are 
called nonlocal atoms to reflect the fact that their energy 
depends on the positions of their neighbors in addition 
to their own position. In contrast, the energies of atoms 
experiencing a smooth deformation field on the atomic 
scale are computed based on the deformation gradient 
in their vicinity as befitting a continuum model. These 
atoms are called local atoms because their energy is based 
only on the deformation gradient at the point where it is 
computed. The total energy E to t (which for a classical 
system can be written as E to t = 2i=i ^ii with Ei the 
energy of atom i) is approximated as 
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The total energy has been divided into two parts: an 
atomistic region of iV nl nonlocal atoms and a continuum 
region of iV loc local atoms (iV nl + iV loc = N r ). The calcu- 
lation in the atomistic region is identical to that in fully 
atomistic methods with the energy of the atom depend- 
ing on the coordinates {q} of the surrounding repatoms. 
However, in the coarse-grained continuum region each 
repatom can represent a large region of rii atoms on the 
atomic scale. Rather than depending on the positions 
of neighboring atoms, the energy of a local repatom de- 
pends on the deformation gradients {F} characterizing 
the finite strain around its position. The basic assump- 
tion employed is the Cauchy-Born rule which relates the 
continuum deformation at a point to the motion of the 
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FIG. 1: QCDFT model for an edge dislocation showing the 
three domains in the model (see text). The x, y and z axes are 
along [110], [111] and [112], respectively. The model contains 
about 990 nonlocal atoms of which 84 are in the DFT region, 
3700 local atoms and 9300 finite elements. All atoms are ini- 
tially displaced according to the anisotropic elastic solution of 
the dislocation with the boundaries held fixed to these values 
during the relaxation phase. Periodic boundary conditions 
are applied along the dislocation line (z) direction. 



atoms in the underlying lattice represented by this point. 
To obtain the necessary deformation gradients, a finite el- 
ement mesh is defined with the representative atoms as 
its nodes. It is important to note that the calculations 
of S* oc ({F}) in the continuum regions make use of the 
same interatomic potential used in the nonlocal atom- 
istic region. This makes the passage from the atomistic 
to continuum regions seamless since the same material 
description is used in both. This seamless description 
enables the model to adapt automatically to changing 
circumstances, for example the nucleation of new defects 
or the migration of existing defects. The adaptability 
of QC is one of its main strengths, which is missing in 
many other multiscale methods. A consequence of the 
partitioning into local and nonlocal regions and the exis- 
tence of a well-defined total energy for the entire system is 
the presence of non-physical ghost forces at the interface. 
These can be eliminated by self-consistent application of 
dead load corrections 0- 

The original QC formulation assumes that the total 
energy can be written as a sum over individual atom 
energies. This condition is not satisfied by quantum 
mechanical models. To address this limitation, in the 
present QCDFT approach the material of interest is par- 
titioned into three distinct types of domains (see Fig.^l: 
(1) a nonlocal quantum mechanical DFT region (region 
I); (2) a nonlocal classical region where Embedded- Atom 
Method (EAM) ;Sj potentials are used (region II); and 
(3) a local region that employs the same EAM potentials 
as region II (region III). The total energy of the QCDFT 
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where E[l + II] is the total energy of regions I and II to- 
gether (the assumption here is that region I is embedded 
within region II). The coupling between regions II and 
III is achieved seamlessly via the QC formulation, while 
the coupling between regions I and II is accomplished by 
a scheme recently proposed by Choly et al. Based 
on this coupling strategy, E[l + II] can be written as 

E[I + II] = E DFT [I] + Seam [H] + E int [I, II] , (3) 

where -Edft [I] is the energy of region I in the absence of 
region II computed using the DFT model, Seam [II] is the 
energy of region II in the absence of region I computed 
using the EAM model, and S lnt [I,II] represents a formal 
interaction energy added to give the correct total energy. 
The interaction energy between the two subsystems can 
be rewritten as: 

E' mt [l, II] = E[l + II] - E[l] - E[II], (4) 
= Seam [I + H] - Seam [I] - S E am [II] . 

The first equation serves as a general definition of the in- 
teraction energy whereas the second equation represents 
one particular implementation of S lnt , which is used in 
this work. Eq. 0] is not contradictory to Eq. |3 because 
EAM has its root in DFT and the EAM energy can be 
viewed as an approximation to the DFT energy. Dif- 
ferent combinations of quantum mechanical and classi- 
cal atomistic methods other than DFT/EAM may also 
be implemented || . The great advantage of the present 
implementation is its simplicity. It demands nothing be- 
yond what is required for a DFT calculation and an EAM 
QC calculation. Furthermore, by substituting Eq.0]into 
Eq. |3I we arrive at 

E[I + II] = S DFT [I] - Seam [I] + S E am [I + H] . (5) 
The forces on the EAM atoms in region II are then 

^ DFT dE EAM [Wl]^S^({F}) 
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where qj 1 are the Cartesian coordinates of atom i in re- 
gion II. It is clear from this equation that the forces on 
the atoms in region II are identical to those that would be 
obtained from a fully-classical QC calculation. The same 
applies to the region III atoms, that is, as far as forces 
are concerned, regions II and III behave as though the 
entire model were classical. This is a very desirable prop- 
erty in terms of achieving a seamless coupling between 
region I and the rest of the model. At the same time, 
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the forces on the DFT atoms in region I will have con- 
tributions from both DFT atoms and the nearby EAM 
atoms in region II. The error in forces on the DFT atoms 
due to the coupling is thus given by the difference be- 
tween calculated forces with DFT and EAM on these 
atoms. To minimize this error, we propose to use a class 
of interatomic potentials which are generated by match- 
ing the forces obtained from the EAM method to those 
from DFT calculations flHflll ]. Another important prac- 
tical advantage of the present QCDFT method is that, 
if region I contains many different atomic species while 
region II contains only one atom type, there is no need to 
develop reliable EAM potentials that can describe each 
species and their interactions. This is because if the var- 
ious species of atoms are well within region I, then the 
energy contributions of these atoms are canceled out in 
the total energy calculation (the last two terms in Eq.[5j|. 
This advantage renders the method particularly useful in 
dealing with impurities, which is an exceedingly difficult 
task for empirical potential simulations. 

The equilibrium structure of the system is obtained by 
minimizing the total energy in Eq. [21 with respect to all 
degrees of freedom. Because the time required to evalu- 
ate -Et>ft[I] is considerably more than that required for 
computation of the other EAM terms in E^ DFT , an al- 
ternate relaxation scheme turns out to be rather efficient. 
The total system can be relaxed by using the conjugate 
gradient approach on the DFT atoms alone, while fully 
relaxing the EAM atoms in region II and the displace- 
ment field in region III at each step. Similar to Choly et 
al. 0, an auxiliary energy function can be defined as 
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which allows for the following relaxation scheme: (i) Min- 
imize E® t DF with respect to the atoms in regions II 
({q 11 }) and the atoms in region III ({q 111 }), while holding 
the atoms in region I fixed; (ii) Calculate £'^ DFT [{q}], 
and the forces on the region I atoms; (iii) Perform a step 
of conjugate gradient minimization of E'; (iv) Repeat 
until the system is relaxed. In this manner, the number 
of DFT calculations performed is greatly reduced, albeit 
at the expense of more EAM and local QC calculations. 
A number of tests have shown that the total number of 
DFT energy calculations for the relaxation of an entire 
system is about the same as that required for DFT re- 
laxation of region I alone. Further computational speed- 
up can be achieved for the DFT calculations by using 
converged electronic charge density and wave functions 
from the previous step, so that the charge (potential) 
self-consistency can be reached faster for the next DFT 
calculation because the atomic relaxation is usually very 
small between two consecutive DFT moves. 

In the remainder of the paper, we apply the present 
QCDFT approach to study the core structure of an edge 
dislocation in Al in the absence and presence of H impuri- 
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FIG. 2: Dislocation core structures obtained from the EAM- 
based QC (top) and the present QCDFT method (bottom). 
The black circles are atoms. The contours correspond to out- 
of-plane (z) displacement in A. The contours clearly indicate 
the splitting of the dislocation. The atoms within the black 
box in the bottom panel are DFT atoms. The finite element 
mesh serves no other purpose in this nonlocal atomistic region 
other than as a guide to the eye to help visualize deformation. 



ties. We chose this system as an example because results 
from both experiments and simulations are available for 
comparison. The QCDFT model for an edge dislocation 
with a Burgers vector f [110] (a = 3.97 A) is presented in 
Fig. ^ Convergence tests on the size of region I indicate 
that a DFT box of 30 A x 9 A x 4.86 A (84 DFT atoms) 
is sufficient to capture the dislocation splitting behavior 
accurately; hence the following calculations are all based 
on this DFT box. A force-matching potential for Al 0] 
was used for EAM calculations. The DFT calculations 
were performed by using the plane- wave pseudopoten- 
tial VASP code nl for a cluster with 8 A vacuum m 
both the x and y directions. The energy cutoff for pure 
Al and Al+H is 129 eV and 200 eV, respectively. We 
find that 10 k points along the one-dimensional Brillouin 
zone are adequate for good convergence. Fig. [2] presents 
the simulation results for both a standard EAM-based 
QC calculation and the QCDFT method, showing the 
dissociation of the edge dislocation into two equivalent 
60° Shockley partials. The splitting distance (obtained 
from an analysis of the displacement jump across the slip 
plane) in the standard QC calculation is 15.4 A, whereas 
the splitting distance obtained with QCDFT method is 
5.6 A, a value very close to the experimentally observed 
value of 5.5 A 13]. This result demonstrates that even 
for a simple metal like Al which should be the best candi- 
date for use of an EAM potential, a quantum mechanical 
calculation is necessary to obtain correct results. 

The most important advantage of QCDFT approach, 
however, is that it allows the study of impurity effects 
on mechanical response, an impossible task for simpler 
empirical potentials. Fig. [31 shows the effect of adding 
one column of H impurities at the dislocation core. The 
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FIG. 3: QCDFT dislocation core structure in the presence of 
a column of H impurities. The circles are Al atoms (black) 
and H atom (white). Contour significance is the same as in 
Fig. H The black lines are a guide to the eye, indicating 
atomic planes. 
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FIG. 4: Charge density distribution in region I in the absence 
(top) and in the presence of one (middle) and two H impu- 
rity atoms (bottom). The blue spheres are Al atoms and the 
red spheres are H atoms. The gray iso-surfaces illustrate the 
charge density distribution at 0.28 electrons/A 3 . Electron 
density values range from to 0.30 electrons/A 3 . 



of H atoms, charge accumulation develops at these H 
atoms as the H impurities attract the valence electrons 
from the Al atoms and become negatively charged. The 
covalent bonding across the slip plane between Al atoms 
is disrupted by the H atoms, and at the same time, ionic 
bonding between the oppositely charged H and Al ions is 
developed. The fact that the directional covalent bonds 
are replaced by more homogeneous ionic bonds near the 
core leads to the wider dislocation core seen in Fig. 

In summary, we have introduced a multiscale modeling 
approach which concurrently couples quantum mechan- 
ical, classical atomistic and continuum mechanics sim- 
ulations, in a unified fashion for metals. Our QCDFT 
method provides a useful framework for multiscale mod- 
eling of metallic materials because it does not require 
the existence of localized covalent bonds for computing 
the coupling energy as all other multiscale methods do 
[J Furthermore, this approach is completely 
general and versatile: it can be applied to diverse ma- 
terials problems, such as dislocations, cracks, surfaces, 
and grain boundaries. Finally, the automatic adaption 
feature of the QCDFT method allows the DFT and/or 
EAM region to move and change in response to the cur- 
rent deformation state, when for example, defects are be- 
ing nucleated in an otherwise perfect region. To demon- 
strate the unique strength of this method in dealing with 
impurities, we have applied it to study H-dislocation in- 
teractions in Al. 

This research was partly supported by an award from 
Research Corporation (GL). ET and GL thank the In- 
stitute for Mathematics and its Applications (IMA) for 
hosting them in the fall of 2004 during which time part 
of this work was done. 



presence of the H atoms results in a spreading of the core 
(the splitting distance is now increased to about 13 A). 
This finding is consistent with the fact, confirmed by ear- 
lier DFT calculations 0] , that H can lower the stacking 
fault energy. The fact that the dislocation becomes wider 
may explain the H-enhanced dislocation mobility that is 
believed to lead to H embrittlement phenomena via the 
enhanced local plasticity theory. A similar core struc- 
ture is also found for two columns of H atoms placed at 
the dislocation core. In order to understand the under- 
lying origin of the H-enhanced dislocation mobility, we 
calculate the electron density distribution at the dislo- 
cation core in the absence and presence of H impurities, 
as shown in Fig. 0] . In the absence of the H impurity, 
the electron bonding is stronger and with a distinct co- 
valent character. The bonding is more directional above 
the slip plane, and it becomes more spherical below the 
slip plane where there are two extra atomic planes, corre- 
sponding to the two partial dislocations. In the presence 
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